import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'simhei'
#读入数据文本
data = np.loadtxt('data.txt')
x = np.zeros((24))
y = np.zeros((24))
#分离X和Y数据
for i in range(24):
    x[i] = data[i][0]
    y[i] = data[i][1]
#输出原始数据散点图
plt.scatter(x,y,label = '原始数据散点图')
plt.legend()
#第一问 选取线性函数y=a0+a1x
a0 = 0
a1 = 0
#基函数为1和x，基函数记为b,内积为b0b0、b0b1等
b0b0 = 0
b0b1 = 0
b1b1 = 0
b0f = 0
b1f = 0
for i in range(24):
    b0b0 += x[i]**0
    b0b1 += x[i]**1
    b1b1 += x[i]**2
    b0f += x[i]**0*y[i]
    b1f += x[i]**1*y[i]
#运用克莱姆法则解二维矩阵
D = b0b0*b1b1 - b0b1*b0b1
D0 = b0f*b1b1 - b0b1*b1f
D1 = b0b0*b1f - b0f*b0b1
a0 = D0/D
a1 = D1/D
#绘制线性图
y0 = a0 + a1 * 1
y1 = a0 + a1 *24
plt.plot([1,24],[y0,y1],label = '线性拟合y=%.3f%+.3fx'%(a0,a1))
plt.legend()
#第二问 抛物拟合 函数设为y=a0+a1x+a2x2
a0 = 0
a1 = 0
a2 = 0
#基函数为1、x、x2，基函数记为b,内积为b0b0、b0b1等
b0b0 = 0
b0b1 = 0
b0b2 = 0
b1b1 = 0
b1b2 = 0
b2b2 = 0
b0f = 0
b1f = 0
b2f = 0
for i in range(24):
    b0b0 += x[i]**0
    b0b1 += x[i]**1
    b0b2 += x[i]**2
    b1b1 += x[i]**2
    b1b2 += x[i]**3
    b2b2 += x[i]**4
    b0f += x[i]**0*y[i]
    b1f += x[i]**1*y[i]
    b2f += x[i]**2*y[i]
#运用克莱姆法则解三维矩阵，采用numpy库函数计算行列式
#构造行列式
matD = np.mat([[b0b0,b0b1,b0b2],[b0b1,b1b1,b1b2],[b0b2,b1b2,b2b2]])
#将b0f,b1f,b2f替换矩阵D各列，得到D0、D1、D2
matD0 = np.mat([[b0f,b0b1,b0b2],[b1f,b1b1,b1b2],[b2f,b1b2,b2b2]])
matD1 = np.mat([[b0b0,b0f,b0b2],[b0b1,b1f,b1b2],[b0b2,b2f,b2b2]])
matD2 = np.mat([[b0b0,b0b1,b0f],[b0b1,b1b1,b1f],[b0b2,b1b2,b2f]])
#使用numpy库函数求行列式
D = np.linalg.det(matD)
D0 = np.linalg.det(matD0)
D1 = np.linalg.det(matD1)
D2 = np.linalg.det(matD2)
#克莱姆法则求未知数
a0 = D0/D
a1 = D1/D
a2 = D2/D
#绘制抛物线
xx = np.linspace(0,24,100)
yy = a0 + a1 * xx + a2 * xx**2
plt.plot(xx,yy,label = '抛物线拟合y=%.3f%+.3fx%+.3fx**2'%(a0,a1,a2))
plt.legend()
#使用温度曲线，可设y=a0+a1cos(w1*x)+a2sin(w2*x)
a0 = 0
a1 = 0
a2 = 0
w1 = 0
w2 = 0
#选取线性空间为1,cosw1x,sinw2x,基函数记为b,内积为b0b0、b0b1等,需要多次求解固写为函数
#w1,w2从可能值遍历，范围取周期附近，2pi/36~2pi/12,步长取（2pi/12-2pi/36）/24
olderror = 9999
curra0 = 0
curra1 = 0
curra2 = 0
currw1 = 0
currw2 = 0
currindex = 0
for i in range(24):
    w1 = 2*np.pi/36+float(i)*np.pi/216
    for j in range(24):
        w2 = 2*np.pi/36+float(j)*np.pi/216
        a0 = 0
        a1 = 0
        a2 = 0
        b0b0 = 0
        b0b1 = 0
        b0b2 = 0
        b1b1 = 0
        b1b2 = 0
        b2b2 = 0
        b0f = 0
        b1f = 0
        b2f = 0
        for k in range(24):
            b0b0 += 1
            b0b1 += np.cos(w1*x[k])
            b0b2 += np.sin(w2*x[k])
            b1b1 += np.cos(w1*x[k])**2
            b1b2 += np.sin(w2*x[k])*np.cos(w1*x[k])
            b2b2 += np.sin(w2*x[k])**2
            b0f += 1*y[k]
            b1f += np.cos(w1*x[k])*y[k]
            b2f += np.sin(w2*x[k])*y[k]
        #运用克莱姆法则解三维矩阵，采用numpy库函数计算行列式
        #构造行列式
        matD = np.mat([[b0b0,b0b1,b0b2],[b0b1,b1b1,b1b2],[b0b2,b1b2,b2b2]])
        #将b0f,b1f,b2f替换矩阵D各列，得到D0、D1、D2
        matD0 = np.mat([[b0f,b0b1,b0b2],[b1f,b1b1,b1b2],[b2f,b1b2,b2b2]])
        matD1 = np.mat([[b0b0,b0f,b0b2],[b0b1,b1f,b1b2],[b0b2,b2f,b2b2]])
        matD2 = np.mat([[b0b0,b0b1,b0f],[b0b1,b1b1,b1f],[b0b2,b1b2,b2f]])
        #使用numpy库函数求行列式
        D = np.linalg.det(matD)
        D0 = np.linalg.det(matD0)
        D1 = np.linalg.det(matD1)
        D2 = np.linalg.det(matD2)
        #克莱姆法则求未知数
        a0 = D0/D
        a1 = D1/D
        a2 = D2/D
        #求平方误差
        error = 0
        for ii in range(24):
            error += (a0*1+a1*np.cos(w1*x[ii])+a2*np.sin(w2*x[ii])-y[ii])**2   
        print('序号',i*24+j,'残差',error)
        if error<olderror:
            olderror = error
            curra0 = a0
            curra1 = a1
            curra2 = a2
            currw1 = w1
            currw2 = w2
            currindex = i*24+j

xxx = np.linspace(0,24,100)

yyy = curra0*1+curra1*np.cos(currw1*xxx)+curra2*np.sin(currw2*xxx)
print(currindex,curra0,curra1,curra2,currw1,currw2)
print(olderror)
plt.plot(xxx,yyy,label = '温度曲线拟合y=%.3f%+.3fcos(%.3fx)%+.3fsin(%.3fx)'%(curra0,curra1,currw1,curra2,currw2))
plt.legend()
plt.savefig('拟合图像.png')
plt.show()